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A slab based long range correction approach for multi-site Lennard-Jones models is presented 
for systems with a planar film geometry that is based on the work by Janecek, J. Phys. Chem. 
B 110: 6264 (2006). It is efficient because it relies on a center-of-mass cutoff scheme and scales 
in terms of numerics almost perfectly with the molecule number. For validation, a series of 
simulations with the two-center Lennard-Jones model fluid, carbon dioxide and cyclohexane 
is carried out. The results of the present approach, a site-based long range correction and 
simulations without any long range correction are compared with respect to the saturated 
liquid density and the surface tension. The present simulation results exhibit only a weak 
dependence on the cutoff radius, indicating a high accuracy of the implemented long range 
correction. 
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1. Introduction 


One of the most important properties of vapor-liquid equilibria that can be deter¬ 
mined by molecular simulation, is the surface tension m . Usually, the properties 
of interfaces are directly sampled in a simulation volume containing both the vapor 
and the liquid phase, separated by an interface. Indirect methods like Grand Equi¬ 
librium j^, NpT plus test particle @] or Gibbs ensemble provide access to the 
bulk properties along the saturation curve in a numerically more efficient manner 
but do not consider interfaces. 

Intermolecular interactions are usually evaluated explicitly up to a specified cut¬ 
off radius, beyond which the interactions are covered by a mean field approach, 
i.e. a long range correction (LRG) which compensates for the cutoff [8l4lCll|. In ho¬ 
mogeneous simulations, the LRG typically only considers the energy and the virial 
III, ll^, while in inhomogeneous systems also the force has to be corrected appro¬ 
priately [i3j[l3]- If a small cutoff radius is used without a LRG, the surface tension 
and other thermodynamic properties are known to deviate significantly from the 
correct values [H, [l^ . 
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For homogeneous systems, typical correction strategies are straightforward, mak¬ 
ing the approximation that the pair correlation function is unity beyond the cutoff 
radius. They may rely on a site-site correction 17| or on center-of-mass correction 
approaches, employing angle averaging [l^ or the reaction field method 
For inhomogeneous configurations, fast multipole methods [H, 22], slab based LRC 
[i 3,[23-[2^ or more complex Ewald summation techniques are used [1^, [ 2 ^. Recent 
implementations of the slab based LRC and the Ewald summation technique yield 


very similar results for planar interfaces IJ, l27| . 


In addition to the LRC approach, the cutoff scheme plays an important role. Eor 
molecular models consisting of several interaction sites, the computational effort 
is much smaller for a center-of-mass cutoff compared to a site-site cutoff. This 
advantage rises with the number of sites. However, the LRC has to be consistent 
with the chosen scheme 1^, 12^ . A site-site cutoff scheme consumes a much larger 
amount of computing time, because every site-site distance has to be evaluated and 
compared to the cutoff radius. E.g., for a pair of carbon dioxide models consisting 
of three Lennard-Jones sites each, the site-site cutoff scheme requires the execution 
of nine distance calculations and if statements for the Lennard-Jones interactions 
during the neighborhood search, while a center-of-mass cutoff requires only one. 
Hence, a center-of-mass cutoff scheme should be preferred. 

In the present work, we combine the slab based LRC approach for inhomogeneous 
systems by Janecek [I^ with the center-of-mass cutoff method by Lustig [I^ . 
which is based on angle averaging, and apply it to molecular models containing 
several Lennard-Jones sites. This combined correction approach is validated for 
planar interfaces with a two-center Lennard-Jones model fluid and two fluid models 
representing carbon dioxide and cyclohexane. 


2. Theory 

The intermolecular pair potential u is usually evaluated in molecular simulation 
explicitly only up to a specified cutoff radius Vc- To correct for the error made by 
this approximation, a LRC has to be applied. The potential energy of molecule i is 
thus separated into the explicitly computed contribution and the LRC contribution 

U,= Y. ( 1 ) 

rij<rc 

Eor systems with planar symmetry, such as a planar liquid film surrounded by 
vapor, it is sufficient to compute the LRC in terms of the coordinate normal to 
the interface, employing a slab-based approach [13, HI. In the present work, this 
corresponds to the y direction. The correction term is then a sum over all 

Ns slabs with respect to the interactions between the molecule i and the 

molecules in slab k 


N, 

k 


( 2 ) 


According to Janecek 
volume 


14|. the correction term is an integral over the slab 


Au 


LRC 

i,k 


2'Kp{yk)Ay / dru(r)r, 


( 3 ) 
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where p{yk) denotes the mean density in slab k and Ay is the slab thickness, cf. 
Figure [H As usual, it was assumed by Janecek 0 that for the radial distribution 
function within a single slab g{r) ~ 1 holds beyond the cutoff radius. 

According to Siperstein, the lower bound of the integration r' has to be selected 
appropriately [ 2 ^ as shown in Figured) If the distance ^ = \yi — yk\ between the 
molecule i and the slab k is smaller than the cutoff radius, the cutoff radius has to 
be used as the lower integration bound, otherwise it is ^ [M, l29|, i.e. 


r'=F' “ (4) 

Itc, instead. 

This definition of r' has to be employed in Eq. ([3]) as well as the analogous expres¬ 
sions for the force and the virial. 



Figure 1. Relevant distances for the LRC approach by Janecek El . If the distance between the molecule 
i and the slab k is smaller than the cutoff radius, the cutoff radius has to be used as the lower integration 
bound, cf. Eqs. m and 0. 


Janecek’s approach yields results that are hardly dependent on the cutoff radius 
for the single-site Lennard-Jones fluid down to Tc = 2.5 a 0,0. It is also suitable 
for multi-site models if the molecular simulation code is based on a site-site cutoff 
scheme. 

However, for molecules consisting of several Lennard-Jones sites, a center-of-mass 
cutoff scheme is more efficient because only the distances between the centers of 
mass have to be evaluated duriM the neighborhood search. In this case, angle av¬ 
eraging as proposed by Lustig [1^ is required for the LRC, because the orientation 
of the molecules cannot be considered explicitly by the LRC. The present study 
introduces such an approach, applying it to the Lennard-Jones potential 


u = Ae _ ^6^-6] 


( 5 ) 


with the energy parameter e and the size parameter a, where s represents the dis¬ 
tance between the interaction sites, which may deviate from the distance between 
the centers of mass r. Three cases have to be distinguished here, cf. FigureE) For a 
given r, the center-center (CC), center-site (CS) and site-site (SS) distances depend 
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on the mutual orientation of the molecules. The term s thus has to be an average 


over all molecular orientations with the same center-of-mass distance r 181]. 




Figure 2. Illustration of the three different cases discussed here. Sites in the center of mass interact with 
each other as a center-center interaction (top), as opposed to the center-site interaction (middle) and the 
site-site interaction (bottom). The distance of the sites from the center of mass of their molecule is denoted 
by T. The dots indicate the center of mass, while the crosses denote the site positions. 


Center-center case 

In the CC case, i.e., for the interaction between Lennard-Jones sites in the center 
of mass, the distance s is equal to the center-of-mass distance r and no angle 
averaging is required, because = r‘^^. For the CC case, the reader is referred 
to Janecek 0 , who derived correction terms for the potential energy, virial and 
force. The present work generalizes Janecek’s approach such that a center-of-mass 
cutoff scheme can be applied to CS and SS interactions with a similar accuracy. 


Center-site case 

In the CS case, a site is not in the center of mass of its molecule, i.e. it is situated 
at a distance r from the center of mass. The CS case does not exist on its own, 
because CC and SS interactions are also always present in such a scenario. The 
angle-averaged value of has been derived by Lustig [I^ 

^ 2 n ^ (r + - (r - 

4rr(n -|- 1) 


where n = —6 or —3, respectively, for the repulsive or dispersive contributions to 
the Lennard-Jones potential. The correction term for the potential energy is then 
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a combination of Eqs. ([3]) and 


= 27rp{yk)Ay / dr4e[(T^^s — (T°s "jr 
Jr' 


2 vre/>(yfc)Ay 


27rep{yk)Aya^ 

3r 


dr 




a 


g (r' + r) ® — (r' — r) ® 3 ^ ^ 

15 ^ 


( 7 ) 


The correction term for the force is obtained in a similar manner 


=-27rp(2/fc)Ay^^ dr^^r 
_ 27rep{yk)Aya‘^C 


rr 


iQ (r' + r) - (r' - r) 4 (r' + r) - (r' - r) 

cr-(T - 


-4J 


( 8 ) 


The correction term for the virial is separated into its normal and tangential con¬ 
tribution. The normal contribution corresponds to the y direction here that is 
perpendicular to the interface, and the tangential contribution corresponds to the 
X and z directions. The term for the virial in normal direction is analogous to the 
force 


AU^^^, = -7rp{yk)Ay dr^^-r 


Tiep{yk)Aya^^ 


2(^2 r 


rr 


10 d-r) — (r' — r) 4 (r' -|- r) ^ — (r' — r) 

a - a - 


-4 


(9) 


The term for the tangential virial is slightly more complicated 

du 


An^Rc = _ 


1 

-jT^p{yk)d^y J dr 


dr r 


7rep(i/fc)Ai/cr^ 

2 rr' 

'^ep{yk)Aya^ 

3r 


10 + t) - (r' - r) 4 (r' + r) ^ - {r' - t)' 

a --- a --- 


(r^ - 


gir' + t) ® — (r' — r) ® 3 ^ “ t) 


-3 


15 


( 10 ) 


Site-site case 

In the SS case, the correction terms are of similar form. Both sites are not in the 
center of mass of their molecule, i.e. they are separated from it by the distances 
Ti and T 2 , respectively. The corresponding expression for has also been derived 
by Lustig 0 


2n — {r + (r — r+) 

8 rrir 2 (n-|-l)( 2 n-|-3) 


2nH-3 


( 11 ) 
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with T+ = Ti + T 2 and r_ = ri — r 2 . The correction terms for the potential energy, 
virial and force are calculated in the same way as for the CS case 


A LRC _ 7rep(yfc)Aya 


12rir2 


^8(^^ + '^+) ^ — W + T_) ^ —(r' —T_) ® + (r'— r+)' 


30 


-fj^[(r' + r+) ^-(r' + T_) ^ + (/- r+) 


( 12 ) 


LRC _ 7rep(yfc)A?/cj3^ 


Mrr = 


2>TlT2r' 


a 


i(r' + r+) ® —(r' + T_) ® — (r'— r_) ® + (r'— r+) ^ 

15 


i(r' + r+) ^-(r' + r_) ^ - (r'- r_) ^ + (r'- r+) ^ 


— a 


(13) 


An 


LRC _ vre/9(yfc)Ay(T3^2 


6 rir 2 r' 


a 


i(r' + r+) ® —(r' + T_) ® —(r' —r_) ® + (/— r+) ® 

15 


(r' + r+) ^ —(r' + T_) ^ —(r' —r_) ^ + (r'— r+) ^ 


— fj 


(14) 


LRC _ vrep(yfc)Ay(T3 


Anl^^h, = 


12rir2r' 


CJ 


(r'+ r+) ® —(r' + T) ® — (/— r) ® + (r'— r+) 


1-9 


15 


— fj 


(r' + r+) ^-(r' + T) ^ - (r'- r) ^ + (r'- r+) 


-3 


A „.LRC 

(15) 


The corrections for the normal and tangential virial are used for the pressure cal¬ 
culation. The surface tension 7 can be obtained from the difference between the 
normal and tangential contributions to the virial n^r — Hy, which is equivalent to 
the integral over the differential pressure PN — PT 


7 = ^ (Hat - Ht) = / dy {pN - Pt) , 


(16) 


where 2A denotes the surface area of the two dividing surfaces [IJ, [31| . 


3. Simulations 


The above correction terms were implemented in the Isl MarDyn molecular dy- 


33l | for an assessment of the present combination of the methods 
18l | . The equations of motion were solved by a leapfrog 


namics code 

by Janecek |14il and Lustig 
integrator 341 with a reduced time step of At = 0.001 a^/rnje for the two-center 
Lennard-Jones model fluid and a time step of At = 1 fs for the real fluids carbon 
dioxide and cyclohexane. Simulations were conducted in the canonical ensemble 
with A = 16 000 molecules. The liquid phase was in the center of the simulation 
volume surrounded by vapor phases on both sides. The elongation of the simulation 
volume normal to the interface was between 60 and 80 a to limit the influence of 
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finite size effects which may be significant for thin liquid films [30|. A thickness of 
the LRC slabs of Ay ~ 0.1 u was used throughout. The spatial extension of the 
simulation volume in the other directions was at least 20 a to account for capil¬ 
lary waves 3^ 371. For the scaling tests, the length of the simulation volume in y 
direction was varied. The equilibration was conducted for 200 000 time steps and 
the production runs for 800 000 time steps. The statistical errors were estimated 
to be three times the standard deviation of four block averages, each over 200 000 
time steps. 

The employed simulation program I si MarDyn was designed for massively 
parallel high performance computing with systems containing a large number of 
molecules [38|. Accordingly, the implementation of the present LRC approach was 
designed to consume only a small amount of computing time and to scale well with 
the molecule number N as well as with the number of processing units. The scaling 
behavior with respect to the number of processing units was discussed in previous 


publications 

here. 


so that no weak or strong scaling experiments are shown 


4. Results 


A series of simulations for the two-center Lennard-Jones model fluid with an elon¬ 
gation of L = a was carried out for different temperatures. The results are com¬ 


pared with those by Stoll et al. [40(] , who employed the indirect Grand Equilibrium 


method where interfaces are absent, with a cutoff radius of Tc = 5 a. In addition, 
simulation results without any LRC are included here, representing the extreme 
case where long range interactions are completely neglected. To exemplify the ne¬ 
cessity of angle averaging, Janecek’s original approach was applied together with 
the center-of-mass cutoff scheme, although it was designed for a site-site cutoff 
scheme (^ . This is termed site-based approach in the following. 



Figure 3. Density p over y coordinate for the two-center Lennard-Jones model fluid for T = 0.979 c/Zcb 
and Tc = 2.5 cr. Comparison between simulations without LRC (dashed line), the site-based approach 
(dash-dotted line), the present approach (solid line) and the reference values by Stoll et al. [4Q| (dotted 
line). 


Figure [3] shows the density over the y coordinate close to the triple point. The 
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density profile is needed on the one hand for the LRC, on the other hand for 
the calculation of the saturated liquid density, which is compared with results 
from indirect vapor-liquid equilibrium simulations. The density with the present 
approach matches the saturated liquid density from the homogeneous simulations 
by Stoll et al. 
data. 


4fll | , while the other approaches exhibit deviations from the reference 



Figure 4. Saturated liquid density over the cutoff radius for the two-center Lennard-Jones model fluid. 
Comparison between simulations without LRC (diamonds ), t he site-based approach (squares), the present 
approach (circles) and the reference values by Stoll et ah [401] (dashed lines). 

Figure 0 ] shows the saturated liquid density over the cutoff radius for different 
temperatures from near the triple point up to 0.96 Tc, where Tc is the critical 
temperature. The results for the saturated liquid density that were determined 
with the present LRC approach hardly show any dependence on the cutoff radius 
for the lower two temperatures, while the site-based approach and the simulations 
without LRC show significant deviations from the reference saturated liquid den¬ 
sity. Simulations without LRC were only performed for comparison near the triple 
point. At the highest temperature, both LRC approaches exhibit deviations from 
the reference case for small cutoff radii. However, it should be noted that a stable 
liquid film at a temperature of 0.96 Tc is quite challenging to simulate due to the 
divergence of the correlation length at the critical point. 

The two-center Lennard-Jones model fluid with L = cr is a difficult case, because 
of its anisotropy. Figure [5] shows the relative deviations from the reference data 
by Stoll et al. for seven model fluids with a varied elongation L at a low 
temperature close to their triple point. These simulations were carried out with a 
constant cutoff radius of Vc = 2.5 a. As expected, the deviations of the site-based 
approach rise with the elongation of the molecules, but they are much smaller than 
in case of the simulations without LRC. The deviations in terms of the saturated 
liquid density reach 3 % for the site-based approach at the largest elongation L 
= a. For small elongations, as expected, the site-based approach and the present 
approach converge. 

Two multi-center Lennard-Jones fluids were also studied to compare the perfor¬ 
mance of the LRC terms: carbon dioxide (CO 2 ), which was described b y a rigid 
three-site Lennard-Jones model with one superimposed point quadrupole |^, and 
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L ! O 

Figure 5. Relative deviations of the saturated liquid density from the reference values by Stoll et ah 
|4Q|| . Comparison between simulations without LRC (diamonds), the site-based approach (squares) and 
the present approach (circles). All simulations were carried out with a cutoff radius of Vc = 2.5 u and a 
temperature close to the triple point, i.e. T 0.56 Tc. 


cyclohexane (C 6 H 12 ), which was described by a rigid six-site Lennard-Jones model 
42l ]. For carbon dioxide, temperatures from 220 to 280 K were considered, i.e. 
almost from the triple point up to approximately 0.92 Tc. Carbon dioxide was cho¬ 
sen because it is similar to the two-center Lennard-Jones model fluid. Moreover, 
it is a combination of the CC, CS and the SS cases. The point quadrupole was 
assumed to haye no preferred orientation beyond the cutoff radius, which yields a 
vanishing LRC contribution. Merker et al. ^ used an indirect simulation method 
without the presence of interfaces and a cutoff radius of at least 7.1 cr in terms of 
the Lennard-Jones parameter a of the oxygen atoms. Figure [6] shows the results 
for the saturated liquid density. The results are similar to the two-center Lennard- 
Jones model fluid, i.e. the saturated liquid density is almost independent on the 
cutoff radius with the present approach. 

For cyclohexane, three different temperatures were studied for a comparison be¬ 
tween the results with different LRC approaches and the reference data by Merker 
et al. [ 4 ^. Figure [7] shows the results for 330, 415 and 500 K. Merker et al. 42| 


also used an indirect simulation method without the presence of interfaces and a 
cutoff radius of at least 4.3 a. 

Because cyclohexane is a much larger molecule than the others considered in this 
work, where all sites have a distance of approximately 0.52 to 0.54 a from the center 
of mass, it is obvious that a cutoff radius of 2.5 a is insufficient. Nonetheless, even 
simulations with a cutoff radius of 3 cr yield good results in terms of the saturated 
liquid density. Only for a temperature of about 0.9 the cutoff radius must be 
larger. 

Another important property of vapor-liquid equilibria is the surface tension. 
For the lowest temperature of the fluids discussed above, the surface tension was 
determined with the three different approaches. Figure [8] shows the surface tension 
over the cutoff radius for the two-center Lennard-Jones model fluid, carbon dioxide 
and cylcohexane. The number of time steps was enlarged to four million to reduce 
the statistical uncertainties and better identify systematic deviations. 

The dependence of the surface tension on the cutoff radius is similar to the depen- 
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Figure 6. Saturated liquid density over the cutoff radius for carbon dioxide. Comparison between simula¬ 
tions without LRC (diamonds), the site-based approach (squares), the present approach (circles) and the 
reference values by Merker et al. (dashed lines). 



Figure 7. Saturated liquid density over the cutoff radius for the cyclohexane. Comparison between simu¬ 
lations without LRC (diamonds), the site-base approach (squares), the present approach (circles) and the 
reference values by Merker et al. (dashed lines). 


dence of the density on the cutoff radius. The present approach shows hardly any 
influence of rc on the surface tension, as opposed to the other discussed approaches, 
which exhibit a significant cutoff radius dependence. 

Furthermore, a simulation series with a single processing unit (Intel Xeon E5- 
2670) with a varying number of two-center Lennard-Jones molecules with an elon¬ 
gation L = a was carried out. The chosen temperature T = 0.979 e/ks is close 
to the triple point of this fluid [d^. Figure [9] shows the computing time for 100 
time steps in the canonical ensemble. Due to the underlying linked-cell algorithm 


[43l . l44l | , the computing time for the explicitly evaluated interactions scales almost 
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Figure 8. Surface tension over the cutoff radius for the two-center Lennard-Jones model fluid, carbon diox¬ 
ide and cyclohexane. Comparison between simulations without LRC (diamonds), the site-based approach 
(squares) and the present approach (circles). The temperature was T = 0.979 e/k-Q for the two-center 
Lennard-Jones model fluid (top), 220 K for carbon dioxide (center) and 330 K for cyclohexane (bottom). 


perfectly with the molecule number N. Only for small systems below N ~ 10'^, 
the LRC does not perfectly scale with the molecule number, because the number 
of slabs does not correlate with it, but rather with the length of the simulation 
volume in y direction. However, even in this case, the computational effort for the 
LRC is more than one order of magnitude smaller than for the explicitly evaluated 
interactions. 
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N 

Figure 9. Computing time for 100 time steps with a single processing unit over the molecule number N. 
The circles correspond to the computing time for the explicitly evaluated interactions. The computing time 
for the LRC (squares) is considerably smaller. 

5. Conclusion 


In this work, a new slab based LRC approach for inhomogeneous systems with 
planar interfaces was presented. It was applied to molecular models consisting of 
several Lennard-Jones interaction sites, employing a center-of-mass cutoff. The 
center-of-mass cutoff scheme is numerically more efficient than a site-site cutoff 
scheme, but it requires a more demanding LRC. The LRC by Janecek that 
is based on the site-site cutoff scheme was generalized to the center-of-mass cutoff 
scheme with the angle averaging method by Lustig (I^ . The influence of the LRC 
on the saturated liquid density and the surface tension was studied. The present 
LRC approach yields very good results for both properties and shows only a weak 
dependence on the cutoff radius. It is numerically efficient, consumes only a small 
amount of computing time and scales well for systems with very large numbers of 
molecules. 
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